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ABSTRACT 

In the present paper, we give the formalism for further developing a fully three-dimensional HZETRN code using 
marching procedures but also development of a new Green’s function code is discussed. The final Green’s function 
code is capable of not only validation in the space environment but also in ground based laboratories with directed 
beams of ions of specific energy and characterized with detailed diagnostic particle spectrometer devices. Special 
emphasis is given to verification of the computational procedures and validation of the resultant computational 
model using laboratory and spaceflight measurements. Due to historical requirements, two parallel development 
paths for computational model implementation using marching procedures and Green’s function techniques are 
followed. A new version of the HZETRN code capable of simulating HZE ions with either laboratory or space 
boundary conditions is under development. Validation of computational models at this time is particularly 
important for President Bush’s Initiative to develop infrastructure for human exploration with first target 
demonstration of the Crew Exploration Vehicle (CEV) in low Earth orbit in 2008. 


INTRODUCTION 

Improved spacecraft shield design requires early entry of radiation constraints into the design process to 
maximize performance and minimize costs. As a result, we have been investigating high-speed computational 
procedures to allow shield analysis starting with preliminary design concepts on to high fidelity final design models 
(Wilson et al. 2003). Of particular importance is the need to implement probabilistic models to account for design 
uncertainties (Wilson et al 2004a) in the context of optimal design processes (Qualls et al. 2003). In this report, we 
will discuss the progress towards a full three-dimensional and computationally efficient deterministic transport code 
for which the current HZETRN evaluates the lowest order asymptotic term in the context of a numerical marching 
procedure. Transport code verification was accomplished through convergence testing and benchmarking using 
restricted solutions with methods of characteristics (Wilson et al. 1988a) and also with Heaviside operational 
transform methods (Ganapol et al. 1991). The marching procedure is adequate for boundary conditions found in the 
space environment allowing space flight validation for nearly two decades using the Space Transportation System 
(STS) and the International Space Station (ISS) measured data. HZETRN is the first deterministic solution to the 
Boltzmann equation allowing field mapping within the ISS in tens of minutes using standard Finite Element 
Method (FEM) geometry common to engineering design practice (VerHage et al. 2002). This coupling enables 
development of integrated multidisciplinary design optimization methods. 
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In the present paper, we give the formalism for further developing a fully three-dimensional HZETRN code 
using marching procedures. In addition we describe the development of a new code based on Green’s functions 
which is capable of not only validation in the space environment but also in the laboratory with directed beams of 
ions of specific energy and characterized with detailed diagnostic particle spectrometer devices (Tweed et al. 2005, 
Walker et al. 2005). A new version of the HZETRN code capable of simulating HZE ions with either laboratory or 
space boundary conditions for validation is under development. The computational model consists of combinations 
of physical perturbation expansions based on the scales of atomic interaction, multiple nuclear scattering, and 
nuclear reactive processes with use of asymptotic/Neumann expansions with non-perturbative corrections. This 
new code contains energy loss with straggling, nuclear attenuation, nuclear fragmentation with energy dispersion 
and downshifts. Off-axis dispersion with multiple scattering and light ion and neutron 3D propagation with full 
angular dependence are under preparation. A benchmark for a broad directed beam for 1 A GeV iron ion beams 
with 2 A MeV width and four targets composed of polyethylene, polymethyl metachrylate, aluminum, and lead of 
varying thickness from 5 to 30 g/cm 2 are also given. The benchmark quantities will be dose, track averaged linear 
energy transfer (LET), dose averaged LET, fraction of iron ion remaining, and fragment energy spectra after 23 
g/cm 2 of polymethyl metachrylate. Flight and ground-based validation of the current codes and their associated 
database are discussed (Wilson et al. 2005a). 

Our early interest in deterministic HZE transport code development paralleled the experimental studies of 
Walter Schimmerling et al. (1986) at the Lawrence Berkeley Laboratory that involved solving the Boltzmann 
transport equation for mono-energetic ion beams in the context of the continuous slowing down approximation and 
velocity conserving interactions (Wilson et al. 1984). The main computational limitation was found to be the 
inadequacy of available nuclear data ascertained by comparing computational results to the ionization data for a 
broad beam of 20 Ne ions (Wilson et al. 1984, Shavers et al. 1993). Emphasis was soon overtaken by the need to 
establish the scope of the galactic cosmic ray (GCR) protection problem, and marching procedures were used to get 
first order estimates of shielding requirements (Wilson et al. 1986, 1991). Validation of the new computational 
marching model against atmospheric air shower data again pointed to the inadequacy of the available nuclear data 
(Wilson and Badavi 1986) and development of a semi-empirical nuclear model followed, leading to reduced 
computational model errors (Wilson et al. 1987a, 1987b). HZETRN code verification was through numerical 
convergence and benchmarking with analytical solution (Wilson and Badavi 1986) and Heaviside methods 
(Ganapol et al. 1991). The next decade emphasized spaceflight validation (including human phantom testing) of the 
marching solution and nuclear model improvements (Shinn et al. 1998, Cucinotta et al. 1998, Wilson et al. 1995a, 
Badhwar et al. 1992, 1994, 2001, 2002). These codes were configuration controlled and accepted by OSHA and 
NASA STD-3000 (2005) for use in radiation safety and design. 

Advanced solution methods of the Boltzmann equation continued to develop (Wilson et al. 1994a,b) but 
only slowly after NASA interest in deterministic transport code and nuclear model development waned in 1995 in 
favor of Monte Carlo methods (Armstrong and Colburn 2001, Pinsky et al. 2001) and their associated limitations 
(Wilson et al. 2002a). Recent renewed interest within NASA for deterministic code development to service the 
CEV design requirements is giving new emphasis to improved solution methods, but with no new nuclear modeling 
activity (Wilson et al. 2004b). As a result, current developments utilize the older semi-empirical NUCFRG2 model 
(Wilson et al. 1995a) with energy downshifts and momentum dispersion (Tripathi et al. 1994, Tweed et al. 2005). 
Comparison of NUCFRG2 nuclear fragmentation database with other models and iron ion beam experiments are 
given by Zeitlin et al. (1997) and Golovchenko et al (2002). In the current report, we will first review the current 
status of code development with emphasis on future needs. Only formal development will be given here, and the 
referenced material contains the details of the computational methods development. Within current limitations, we 
will present methods of code verification and recent spaceflight validation based on the marching procedures and 
laboratory validation using a newly developed code based on three-term Neumann expansions with nonperturbative 
corrections (Wilson et al. 1994a, 1994b) above the third Neumann term (Tweed et al. 2005). 

DETERMINISTIC CODE DEVELOPMENT 

The relevant transport equations are the linear Boltzmann equations derived on the basis of conservation 
principles (Wilson et al. 1991) for the flux density (f)j(x, Q, E) for particle type j as 

G*V(t)j(x,Q,E)=YS Ojk(£2,n',E,E') ^(x&'p') dQ' dE' - Oj(E) (j)j(x,Q,E) (1) 
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where Oj(E) and Ojk(£2,£2',E,E') are the shield media macroscopic cross sections. The Ojk(£2,£2',E,E') represent 
all those processes by which type k particles moving in direction £2 / with energy E' produce a type j particle in 
direction £2 with energy E (including decay processes). Note that there may be several reactions that produce a 
particular product, and the appropriate cross sections for equation (1) are the inclusive ones. Exclusive processes 
are functions of the particle fields and may be included once the particle fields are known. The total cross section 
<Jj (E) with the medium for each particle type is 


<jj (E) - Cj y at (E) + 0) } el (E) + Ojj (E) 


( 2 ) 



where the first term refers to collision with atomic 
electrons, the second term is for elastic scattering on the 
nucleus, and the third term describes nuclear reactions 
where we have ignored the minor nuclear inelastic 
processes (excited single particle states, except for low 
energy neutron collisions). The corresponding 

differential cross sections are similarly ordered. Many 
atomic collisions (~10 6 ) occur in a centimeter of 
ordinary matter, whereas ~10 3 nuclear coulomb elastic 
collisions occur per centimeter, while nuclear reactive 
collisions are separated by a fraction to many 
centimeters depending on energy and particle type. We 
include in the Oj y r(E) term the nuclear decay processes. 

Solution methods first use physical perturbations based 
on the ordering of the cross sections with the frequent 
atomic interactions as the first physical perturbation 
with special methods used for neutrons for which atomic 
cross sections are zero. The first physical perturbation to be treated is the highly directed atomic collisions with 
mean free paths on the order of micrometers as observed in nuclear emulsion. The usual approximation is the 
continuous slowing down approximation leading to well specified range-energy relations as shown in Fig. 1, but 
neglects the energy straggling that will be included in the present treatment. The next term is the highly directed 
multiple coulomb scattering and is usually neglected in many models but is of great importance in understanding 
the transport of unidirectional ion beams leading to beam divergence and so is treated in detail. The remaining 
nuclear reactive processes have been the main attention in past code developments. We will now consider the 
formal development of the relevant equations for further consideration. 

We rewrite equation (1) in operator notation 
by defining a vector array of field function as 

O = [tyx, £2, E) ] 
the drift operator 


Energy (MeV/amu) 

Fig. 1 Range of ions in aluminum. 


D = [£2*V ] 

and the interaction operator 
I = [IS Gjk(£2,£2\E,E) d£2 ' dE' - Oj(E)] 


(3) 

(4) 

(5) 


with the understanding that / has three parts 
associated with atomic, elastic, and reactive 
processes as given in equation (2). Equation (1) is 
then rewritten as 


[D-I at -I el ]*0= I r *0 


( 6 ) 


where the first two physical perturbation terms are 
shown on the left-hand side and are represented by 
diagonal operators. The first order physical 
perturbation for atomic processes is the solution of 



Fig. 2 FWHM of 49.1 MeV protons. 
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[D-IJ»0= 0 


(7) 


t.: 


Using the moments methods is approximated by 
O (z,E) = exp[- (E - (E (z))f / (2 s (z))]l ( V(2 n) s(z)) (8) 

where the array of mean residual energies (E(z)) is 
calculated from the range-energy relations given in Fig. 1 
and the energy deviation due to straggling s(z) are 
evaluated using a second order Green’s function (Wilson 
et al. 2002b) and related to the full width at half maximum 
(FWHM) as shown in Fig. 2 for 49.1 MeV protons 
penetrating aluminum targets. As a result of energy 
straggling as given in equation (8), the unique relationship 
of penetration depth and particle energy as shown in Fig. 1 
is no longer true. The fluence of penetrating protons in 
aluminum is represented by equation (8) and shown in 
Fig. 3 where the mean range is given by the curves in Fig. 
1. Although the straggling correction for the uncollided 
beam with monoenergetic boundary functions is 
important, it is negligible in the higher order terms 
compared to the energy dispersion in nuclear reaction 
products. 

The second physical perturbation term is the 
multiple coulomb scattering by the atomic nucleus as a 
solution of 


[D-I el ]*O = 0 (9) 

and represented by Rutherford scattering modified by 
screening of the nuclear charge by the orbital electrons 
using the Thomas-Fermi distribution for the 
approximation of atomic orbitals. We will utilize the 
multiple scattering solutions of Fermi given by 

0 (z,r,6 n 0 t ) -[^3 w 2 /27iz 2 ] 2 xexp[-w 2 ( 9 2 /z)] 
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Fig. 3. Proton fluence near end of range in aluminum at three 
energies. 



Fig. 4. Characteristic scattering length of various ions in 
aluminum. 


xexp[ - w 2 ( 6 2 /z-3r 6 r /z 2 +3r 2 /z 3 ) ] (10) 

where z is the longitudinal distance, r is the lateral 
distance, 6 r is the angle to the longitudinal axis in 
the r,z-plane, 9 t is the propagation angle 
perpendicular to the r,z-plane and w 2 the array of 
appropriate angular diffusion coefficients. Strictly 
speaking, the solution applies only over intervals 
for which the variation in ion energy is small. It 
follows that the mean square transverse angle 
(understood as a differential quantity) in traveling 
a distance dz is given from equation (10) as 

(dr)=2dz/w 2 ( 11 ) 

Conversely, one finds for a uniform nuclear 
charge distribution shielded by a Thomas-Fermi 
atomic structure 

<9 r 2 ) = dz (Z P E s /ppc) 2 /X 0 (12) 

so that 

w 2 = 2X 0 (P pc/Z P E s ) 2 (13) 


Multiple Coulomb Scattering of 600 MeV/cmu Fe +IB on Aluminum 



Fig. 5 Multiple Coulomb scattering of 600 A MeV iron ions in a 5 g/cm 2 
aluminum target. 
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where 

E s = V( 4n He /a) = 21.2 MeV (14) 

with a the fine structure constant, p the ion momentum 
array, the array /J of ion speeds relative to the speed of 
light c, Z P the array of projectile charges, and X 0 is the 
electron radiation length in the material (Fig. 4). The 
electron radiation length (g/cm 2 ) is given by 

Xq 1 = 4a (Na/A) Z T 2 r 0 2 ln(181Z T - 1/3 ) (15) 

with Avagadro’s number N A , A the molecular weight of 
the target material, and r 0 the classical electron radius. We 
will be using Schimmerling and coworkers modifications 
to Fermi’s width formula (Schimmerling et al. 1986, 

Wong et al. 1989). The effects of energy variation on the 
variance in angle is given by 

(Or 2 ) = I [ (Z p E s /(3pc) 2 /X 0 ] dE/S/E) (16) 

where the integral is over the energy range, and S/E) is 
the stopping power. Multiple scattering played a critical role in prior experimental validation of the transport 
solutions (Shavers et al. 1990,1993). An example of multiple Coulomb scattering is given in Fig. 5 showing the 
emerging ion angular distribution on the beam axis and off the beam axis for 600 A MeV iron ions in aluminum. 
The iron ions on the beam axis remain highly peaked in the forward direction while those that have scattered off 
axis exhibit an angular divergence related to the longitudinal distance z and lateral distance r as 6 r = 3 r/(2 z). In all 
cases, the angular dispersion is small except at very large distances z and clearly will be mainly important only for 
the uncollided beam ions. To a first approximation, the off axis components are taken as directed along the central 
angle 6 r = 3 r/(2 z) with dispersion width of in-plane and out-of-plane given by equation (16) and treated in an 
asymptotic expansion about the central angle. Note that uncharged particles have no Coulomb scattering, and the 
main contribution to equation (9) is the nuclear elastic scattering and will be discussed later in the development. 

The third order physical perturbation involves the relatively rare nuclear reactive processes represented by 
the operator l r of equation (6). The probability of nuclear reaction in coming to rest is given in Fig. 6. It is seen 
from the figure helium ions have unit probability of nuclear collision above 500 A MeV and other ions will only 
reach unit probability of nuclear reaction above 1 A GeV. Although these nuclear processes are relatively rare, they 
represent large transfers of energy resulting in production of additional particles in the spacecraft interior 
environment which are at times more injurious than the particles of the space environment present at the boundaries 
of the spacecraft. The nuclear absorption term is a diagonal contribution and will be shifted to the left hand side 
with other purely diagonal terms while the integral contributions have both off-diagonal and on-diagonal parts and 
will remain on the right hand side as follows 

[D-i at -i el + <j r ]*0 = {S(j r (n,n',E,E')da'dE}*0 =s><z> (17) 

It is the control of secondary processes represented by the right hand side of equation (17) that is a controlling 
factor in reducing the radiation environment within a spacecraft. This is accomplished by choice of chemical 
composition of construction materials that must simultaneously meet other requirements of the service environment 
including loads, thermal, damage tolerance... (Wilson et al. 2002c, 2004a). The elements with off-diagonal nuclear 
reactive differential cross sections can be written in the following form 

o jkj (Q,Q',E,E') = Oj k>iso (E,E')l4n + a M JQ,Q)E,E) (18) 

where the first term is isotropic and associated with lower energy particles produced including target fragments and 
the second term is highly peaked in the forward direction and is associated mainly with direct quasi-elastic events 
and projectile fragmentation products (Wilson 1977, Wilson et al. 1988a). 

Surprisingly, even nucleon-induced reactions follow this simple form and the isotropic term extends to 
relatively high energies (see Fig. 7). For nucleon induced reactions, the following form has been used in early 
versions of FLUKA (Ranft 1980, Wilson et al. 1988b) as follows 



Fig. 6 Probability of nuclear reaction as a function of ion type 
and energy. 
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(19) 


G lk ,(an',E,E') = v jk (E') C7 k ,(E% k (E,E') x g R (Q*Q',E,A T ) 

where V jk ( E') is mean multiplicity and the Ranft factor used in FLUKA is 

g R (Q*Q ' E,A T ) =N r exp[-ffa R ] 7v2>6>0 (20) 

and constant for larger values of production angle 6 , is normalization constant, and X R given by Ranft as 

X R = (0.12 + 0.00036A t /E) (21) 

although new generalized fits are being derived. This separation in phase space will be further exploited in 
computational procedures. 

The heavy ion reaction cross sections have the general form (Wilson 1977) given by 

a jkM (£2,£2',E,E')= o k ,(E') [S jk (E')f jk (£2,£2',E,E’) + \? jk (E’)f jk (£2,£2',E,E')] (22) 

where superscript P refers to projectile fragments and T to target fragments and are largely recognized by their 
separation in phase space. The heavy ion projectile fragment cross sections are further represented by 

(7 jkS JQ,Q',E,E')= o k ,(E')\r jk (E')N, exp[-2mV(E')(l-£2»£2')/£ tjk ] xexp[ - (E + X jk -E?/2 £ 2 ] N(2n£ jk 2 ) (23) 

where X jk is related to the momentum downshift, e jk is related to the longitudinal momentum width, e t jk is related to 
the transverse momentum width, and N t is the transverse normalizing factor. Since the transverse width is small 
compared to the projectile and fragment energy the transverse function is highly peaked about the forward direction 
(£2^£2' ~ 1). The nucleonic component spectral distribution in HZE collisions for quasi-elastic events is 
approximately given as 

C7 jkJ J£2,Q',E,E')= O k JE')\f ]k (E') N n xexp[-2n, V(E')( !-Q -Q ')/£, Jk ]exp [- 1 E-E'\ K /2w 2 f ] (24) 

where R is given in Shavers et al. (2001) as 

R=L3108+0.883exp(-0.0038E')+0.6933exp(-0.001E / ) (25) 

and the width Wf is given as 

w f = 197.3x(1.125A k ) 1/3 /T k (26) 

where T k is the Fermi radius of the projectile ion k. 

Atomic interactions limit the contributions of charged 
particles in the transport process. For example, the 
protons and alpha particles produced in aluminum below 
100 A MeV contribute to the fluence only within a few 
centimeters of their collision source and the heavier ions 
are even more restricted (see Fig. 1). This is an 
important factor in that the transported secondary 
charged particle flux tends to be small at low energies 
and the role of additional nuclear reactions are likewise 
limited (see Fig. 7). 

The reaction cross section is separated by 
equation (18) into isotropic and forward component for 
which equation (17) may be written as coupled 
equations 



Fig. 7 Isotropic and forward neutron spectra produced by 500 MeV 
protons in aluminum. 


[D - I at -I el + a r ]» 0 for = { j a rf JQ,Q',E,E') dQ' dE'}» 0 fur = S rJor » ® fur (27) 

and 

[D - I a , -I e i + c r ]* Oise = {\ 0/£2,£2',E,E') dQ'dE}' O iso + {\a r J£2,£2',E,E') d£2'dE'}» O fo ^S r » O iso + S r , iso » 0 for (28) 

Equation (27) is used to match the boundary conditions and can be written as a Volterra equation (Wilson 1977, 
Wilson et al. 1975, 1991, 1994a, 1994b) and solved as a Neumann series 

'for / E "F G • HyJ 0r *G “I" G • EiyJ 0r *G • n r j or *G +•.,] O&B 


(29) 
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for which the series can be evaluated directly as described elsewhere (Wilson et al. 1994a). 

The cross term in equation (28) gives rise to an isotropic source of light ions of only modest energies and 
neutrons. The high-energy portion of the isotropic spectra arises from multiple scattering effects and the Fermi 
motion of the struck nucleons within the nucleus. The low-energy isotropic spectra arise from decay processes of 
the struck nucleus. Spectral contributions to the Neumann series depend on the particle range and probability of 
surviving nuclear reactions that establish the functional form of the G matrix in equation (29). The second term of 
the Neumann series is proportional to the probability of nuclear reaction that is limited by the particle range as 
discussed above and shown in Fig. 6. It is clear from Fig. 6 that those nuclear reactions for the charged particles 
below a few hundred A MeV are infrequent providing fast convergence of the Neumann series (Wilson and Lamkin 
1975). For the moment, we will neglect the straggling and multiple-elastic processes to simplify the present 
explanation (these provide only minor corrections to space radiation exposures but important in laboratory testing) 
and examine the remaining reactive terms of equation (27). The corresponding Volterra equation is given (Wilson 
1977) by 

<j)j(x, Q, E) = {S/EJP/EJ <pj(r(Q,x), Q, Ej + if/ElE 'AjP/E ')f~f 4x dE "dQ ' a jkf JQ, Q',E',E") 

X<t> k (x+[R j (E)-R j (EOian',E'')}/S ] (E)P j (E) ’ (30) 

where T is the point on the boundary connected to .v along -Q. E y = Rj 1 [p - d + RJ, p is the projection of x onto Q. 
and d is the projection of T onto Q. Equation (30) results from the first term of the Neumann series solution to 
equation (27). In the past we have expanded the angular integral over Q ' asymptotically and implemented as a 
marching procedure (HZETRN, Wilson and Badavi 1986), as a perturbation expansion (Wilson et al. 1984), and by 
non-perturbative approximation (Wilson et al. 1994a) resulting in three distinct methods to evaluate the first order 
asymptotic terms, all of which have had extensive experimental validation (Shavers et al. 1993, Wilson et al. 1998, 
Shinn et al. 1998). Independent of the method used to evaluate the lowest order term, the first correction term is 
found by replacing the fluence in the integrand of equation (27) by the lowest order asymptotic solution as 

<j)j(x, Q, E) = {S/EJP/EJ <t>j(r(Q,x),n,E) + Zf^dE'AjP/E'^L dE"dQ' a jkJ JQ, Q \E',E") 

X (f> ko (x+[Rj(E) - Rj(E')]Q',Q',E")}/ SJEjPJE) ’ (31) 

where (f)j(x, Q, E) is found as an integral over the 
neighborhood of rays centered on Q using the lowest order 
asymptotic solution (j) ko (x,£2',E") along an adjacent ray 
directed along Q'. Note that the boundary condition reached 
along -Q' enters through the lowest order asymptotic 
approximation and the angular integral correction in 
equation (30) is determined by the homogeneity and angular 
dependence of the space radiation and radius of curvature of 
the bounding material as we have shown long ago (Wilson 
and Khandelwal 1974, Wilson 1977). These are the 
determinant factors of the magnitude of the first order 
asymptotic correction that is anticipated to be very small for 
human rated systems (large radius of curvature) in space 
radiation that is homogeneous and isotropic in most 
applications (Wilson et al. 1991, Wilson et al. 1994b). 

In a region of small radius of curvature the specific 
flux components near the site of evaluation will be missing 
contributions along adjacent rays that do not compensate 
losses along the ray on which the solution is evaluated representing the losses due to leakage. (Note, an asymptotic 
treatment of such small angle dependent phenomena is the only useful approach circumventing large discretization 
errors.) This computational procedure is only a small addition to prior code development and will have little impact 
on computational efficiency. The angular dependence of the integral kernel of equation (31) is controlled by the 
forward reactive cross section <jj kfor ({2,{2',E',E") with its highly peaked structure given by equations (19, 23, 24) 
depending on particle type. The angular dependence of the forward peak of fragmenting Ca ions at 100 and 1,000 
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Fig. 8 Normalized transverse components for 
Ca fragmentation. 


A MeV is shown in Fig. 8 as an example. The low-energy ions with limited range have transverse components on 
the order of 10 degrees reducing to a few degrees at high energies. It should be clear that the added divergence 
added by multiple Coulomb scattering of such fragments (Fig. 5) is negligible to the large angular widths of the 
fragmentation event (Fig. 8) further justifying equations (30 & 31). 

Note that the low energy ions have limited range and will contribute little to the transported flux (see Fig. 
1) or nuclear reactions (see Fig. 6). The higher energy ions, with their much longer pathlengths therefore giving 
more important added contributions to the particle fields, are related to only a very small angle of acceptance (few 
degrees) at the boundary. The form of the kernel leads directly to a Gauss-Hermite expansion and evaluation over 
the angle of production. Although the neutron Neumann series for the forward components converge more slowly 
since their contribution to the neutron flux is not limited by atomic interactions, these higher energy neutrons will 
be adequately evaluated by similar procedures. Higher order asymptotic terms can be evaluated with similar 
iteration of equation (31) if required but all indications are that the first such correction will be small for space 
radiation (Wilson and Lamkin 1975). This leaves the diffuse components of neutrons and light ions produced in the 
collision of the forward components and transported by equation (28) to be resolved (see for example, Clowdsley et 
al. 2000). 

The transport from the low-energy neutron and light-ion isotropic sources in equation (28) dominates the 
solution below about 70 A MeV (see Fig. 7). In this region light-ion transport is completely dominated by the 
atomic interaction terms, and only a very small fraction have nuclear reactions making only minor additional 
contributions to the particle fields. This is especially true for the target fragments that can be solved in closed form 
(Wilson 1977, Cucinotta et al. 1990, 1991). The transport solution for the isotropic ion source terms to the lowest 
order perturbation is given by 

<t> jiiso o(x,Q,E) = Zj,! : hlE'A,P,(E')J l: :i 4K dE"d£2'a jk . s J£2,£2’,E',E") <t> kfor (x+ : [R/E) - R/Ey£2,£2',E')/ S/EJP/E) (32) 

and will give highly accurate solutions to equation (28) since very few of the ions will have reactions (see Fig. 6), 
but could be easily corrected using the HZETRN light-ion propagator applied to the diffusive source terms. Note 
the E / integral effectively sums the ion source terms along direction Q from the boundary to x. Also, the nuclear 
survival terms Pj(E) are all near unity (see Fig. 6 showing 1 - Pj(E)). 

The neutrons have no charge and are 
undergoing, at low-energies, mainly elastic but also 
reactive nuclear processes. Although, equation (28) 
exhibits behavior similar to thermal diffusion, there 
are strong differences between thermal (heat 
conduction) and neutron diffusive processes. Thermal 
diffusion at ordinary temperatures has minor leakage 
through near boundaries since radiative processes are 
proportional to fourth power of temperature (in the 
absence of convection) leaving lateral diffusion across 
the media an important process. In distinction, 
neutron diffusion is dominated by leakage at near 
forward and backward boundaries and lateral diffusion 
plays a minor role. Generally, low-energy neutron 
leakage is a dominant process within 15-20 g/cm 2 of 
the bounding surface in most materials. Since human 
rated systems have shields of large radius of curvature 
and small thickness to radius ratio as determined by 
living and working space requirements, it 
approximates a connected system of flat plates for 
which leakage at forward and backward boundaries 
dominates the neutron transport. In this limit, neutron transport simplifies to a sequence of approximate ID- 
transport problems with leakage at the back and forward boundaries and shows reasonable success in comparison 
with experimental flight data (Clowdsley et al. 2000, 2002). In the present development we will consider a 
convergent series of approximations to gauge accuracy of the transport procedures and allow choices of the most 
practical method. 

The first term for diffusive neutron transport uses the lowest order perturbation similar to equation (32) 

given as 



Fig. 9 Spectral dependence of diffuse isotropic neutron 
source term (eq. 33) and the first collision source term 
(eq. 34). 
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(33) 


v °(x,Q,E) = Zl f r l dx'exp[- o n (E)x'] Se "4 dE "dQ ' a jU JQ, Q',E',E") <j> kf J x - x'Q,Q',E") 


where p-d is the distance along -Q from x to the boundary, x' is the distance from x to the source point along -12, 
and o n (E) is the total neutron cross section (elastic plus nonelastic). Note that equations (32) and (33) transport all 
particles (charged and uncharged) associated with the collisions of the forward component & for found as solution to 
equation (27). The second collision term associated with the diffuse charged particle field given by equation (32) is 
negligible, but additional source terms from the lowest-order diffuse neutron solution given by equation (33) 
provides a strong source of diffuse neutrons. The added transport of these neutrons is given by 


[Q*V 


+ a„(E)]<P n J(x,Q,E)=\ O n (Q,Q',E,E')dQ'dE’(P n 


>(x,n',E') + J<7 (Q,Q',E,E’)dQ'dE'6 J'(x,Q',E') (34) 


where (t> niso 1 C x,Q,E) is the remaining diffuse neutron component. The source term to the far right of equation (34) 
have been solved in exact 3D geometry, and the energy spectrum is much degraded from the source term in 
equation (33). Typical spectra of the sources in equations (33) and (34) are shown in Fig. 9. It is clear from the 
properties of the second collision source term that the diffuse spectra of neutrons from this term are highly 
degraded in energy and the methods developed in the nuclear engineering community for reactor applications are 
fully applicable, and our final attention turns to solution of equation (34). Note that the software prepared for 
equation (34) is also applicable to nuclear reactor shielding issues with appropriate source terms. 

The dominant contribution to the low-energy neutron transport in most materials is elastic scattering from 
the media nuclei. In our shell FEM geometry, equation (34) can be written as 


AM + a n (E)]dJ(x,p, E) = 1 o(Q,Q',E,E')dQ'dE'<P n J(x,p',E') + ja 


(Q, Q ’,E,E’)dQ 'dE'<P niso «(x,{2',E') 


(35) 


where q = Q»e x , q' = Q'*e x , the last term on the right is the source term from the forward scattered particle 
collisions. The elastic scattering term in equation (35) has a unique angular contribution p 0 = for a given 

energy transfer (E'~ E) relating to a unique direction q' under the integral. Standard spherical harmonic expansions 
of the cross sections and flux are made to develop a solution. As example, first order Legendre expansion (PI 
approximation) gives the forward and backward flux as 


4>f,b — 0.5 4>o ± 0.75 (|)i 


(36) 


where cf>i are the Legendre coefficients of </> nis J( x,p,E) with approximate transport equation given as we have 
used in the past (Clowdsley et al. 2000) using multigroup and collocation methods. 

The multigroup transport of the diffuse neutrons in multilayered flat plate geometry with variable front and 
back boundaries in which the ions are treated in 
the lowest order asymptotic approximation is the 
current production code HZETRN. The first 
improvement would be to treat the diffuse ion 
components with perturbations from the neutron 
diffuse component giving a complete flat plate 
code in lowest asymptotic order with multigroup 
neutrons and diffuse ion components. To this, we 
will add the first asymptotic correction. This will 
provide NASA with an interim HZETRN 
production code for engineering design process 
development. In addition to the above, we will add 
mesons and electromagnetic cascades in the near 
future (Blattnig et al. 2004, 2005, Nealy et al. 

2002, Anderson et al. 2004). At each step of future 
improvements, the corresponding codes will be 
integrated into the Space Ionizing Radiation 
Environments and Shielding Tools on the 
Framework for Analysis and Collaborative 
Engineering (SIREST/FACE) environment 
(Singleterry et al. 2001) with environmental models for NASA use (the current production code resides there 
already). 



Fig. 10. Fluence profile with marching procedures that are 
indistinguishable from the Heaviside solutions of Ganapol et al (1991). 
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VERIFICATION 


The solution using marching procedures is dependent on step size and energy grid, and convergence to 1 
percent in the first 200 g/cm 2 of material is typically used for shield effectiveness evaluation (Wilson et al. 1991). 
The solution also depends on the isotopes represented in the calculation. Early versions of the code represented 
only the spectrum of 28 charges (H through Ni) plus neutrons. In this case, the fluence of hydrogen ions was the 
sum of proton, deuteron, and triton fluence combined and left the individual isotopes undefined. The marching 
procedures of this code was benchmarked against a Heaviside solution with a simplified stopping power and 
fragmentation parameter set allowing such techniques (Ganapol et al. 1991). The comparison of the marching 
solution and the Heaviside solution is indistinguishable as shown in Fig 10 demonstrating the high accuracy of the 
HZETRN algorithm. 





Tiis-ue dtp Eh., 



Fig. 11. Dose (left panels) and dose equivalent (right panels) for a rigidity spectrum with rigidity parameter of 100 MY incident 
on 20 g/cm 2 iron shield ahead of 30 cm tissue. Penetrating primary contributions (top panels) and secondary contributions 

(lower panels). 
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Monte Carlo methods were the first to 
develop under the Apollo program for spacecraft 
design. The High Energy Transport Code (HETC) 
developed at Oak Ridge National Laboratory for 
space use under NASA contract had at its core the 
Bertini intranuclear cascade model with a strong 
experimental validation effort and the low energy 
neutron data files of the nuclear energy program 
(Bertini et al. 1972, Alsmiller 1967). Although an 
accurate procedure, the computational burden 
even in simple geometry limited their usefulness 
in practical applications. Simpler and less 
accurate computational models were finally used 
in Apollo design and operations as the 
computational burden made HETC impractical. 
Only in the Shuttle operations era were improved 
computational procedures developed allowing 
computational efficiency based on the marching 
procedures discussed herein (Wilson et al. 1988c, 
1991). The HETC code in simple geometry 
allowed verification of these new computational 
procedures and database were reasonably correct. 
One such comparison consisted of the benchmark 
using a solar particle event rigidity spectra based 
with rigidity parameter of 100 MV. The 
comparison of the marching procedure with the 
HETC computation for a 20 g/cm 2 iron shield of 
30 cm of tissue is shown in Fig. 11. It is clear that 
the marching procedures and preliminary database 
give a fair representation of such integral 
quantities as dose and dose equivalent for typical 
solar particle events. 

A benchmark for analytic solution again 
confirmed convergence to better than 1 percent 
(Wilson et al. 1988a). A convergence test by 
adding isotopes representing masses from 1 to 58 
with neutrons indicated convergence to 1 percent 
in dose equivalent for the full galactic cosmic ray 
spectrum and was the standard used until recently. 
Isotope representation for particle fluence 
measurements for iron beams required a 122 
isotope representation (Kim et al. 1994). Recently 
the galactic cosmic ray isotope distribution was 
increased to 170 isotopes and the corresponding 
isotope distribution in the transport model 
includes all such isotopes with expected 
convergence in dose equivalent to be much less 
than 1 percent (Cucinotta et al. 2003). Many more 
years of experience has been gained in validation 
with both ion beams and in spaceflight. By far the 
most important validation of the transport codes 
has been the laboratory validation process where 
both atomic and nuclear interaction details of 
movement through the material are gained. But 



LET, keV/jim 

Fig. 12. Let spectra of a 557 A MeV iron beam measured in nuclear 
emulsion in a layered shield arrangement of 2.24 g/cm 2 of lead followed 
by 2 g/cm 2 aluminum (Shinn et al. 1994). 



Fig. 13. Energy deposition in water by 20 Ne ions at 670 A MeV as 
evaluated using Green’s functions and the Silbergberg et al. (1977) 
nuclear model and the VR correct model. 
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Fig. 14. Integral fluence of fragments produced by 636 A MeV neon ions 
in water of varying thickness. 
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Fig. 15. Charge-changing cross sections for AZ from -1 to -14 for 1.05 

GeV/nucleon 56 Fe incident on H, C, Al, Cu, and Pb targets. The solid 
lines are predictions from NUCFRG2 database generator. 
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spaceflight validation has the overall ability to test adequacy of shield modeling tools, environmental models, and 
transport models. Such a complex test has the disadvantage of ambiguity of errors that can enter the validation in 
different ways. To test environmental model errors, we have relied more on the better defined shielding model of 


the Shuttle in which high resolution 
measurements have been made over long 
periods of time. 


* Expo-inert 

676A MeV^-frC 674A MgV^Ni-N 672AMgV 16 0-» 0 


ION BEAM VALIDATION 

The current status of computational 
procedures are detailed by Tweed et al. 
(2005) with validation of these procedures 
by Walker et al. (2005). In the past we have 
found most useful the direct comparison of 
energy loss spectra either measured in 
nuclear track detectors (Shinn et al. 1994, 
Golovchenko 2002) or in silicon detectors 
(Wilson et al. 1984, Shavers et al. 1993, 
Wilson et al. 1998). An example of 
validation using nuclear track detectors is 
shown in Fig. 12. In this case, the target 
consists of two layers of materials, lead is 
the first target encountered by the beam and 
is followed by an aluminum target. The 
peak on the far right is the penetrating iron 
ions that have either survived penetrating 
the target or have only neutrons removed in 
interactions. The widths of the peaks are 
mainly a result of the resolution of the 
detectors and not a characteristic of the ion 
beam or interactions with the target. 
Validation with electronic detectors began 
with the laboratory tests of neon ion beams 
by Schimmerling and coworkers using a 
water column and ion chambers shown in 
Fig. 13. Also shown in the figure are model 
results of Wilson et al. (1984) in which the 
nuclear model used was a critical issue and 
replacement of the “state-of-the-art model” 
developed by the Naval Research 
Laboratory (Silberburg et al. 1977) required 
correction and finally replacement with 
software using a new nuclear model named 
NUCFRG (Wilson et al. 1987a, 1987b). The 
curves D (0) , D (1) , and D (2) are successive 
perturbation terms using the Wilson et al 
(1984) corrected cross sections. A further 
comparison of the fragment fluence is 
shown in Fig. 14. Note detector efficiency 
problems for the Be and B fragments 
(Shavers et al. 1993). Blind cross section 
testing by Zeitlin et al. (1988) confirms the 
usefulness of the NUCFRG model (Fig. 15). 



Fig. 16. Attenuation of primary beam particles (dash) and secondaries formed 
by neutron removal (solid). 
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Fig. 17. Secondaries produced by removal of one proton. 
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Fig. 18. Secondaries produced by removal of two protons including alpha 


clusters. 





pepfb, cm 

Fig. 19. Secondaries produced by removal of three protons. 
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Further laboratory testing of LET spectra at several 
depths in water using electronic detectors (Shavers et al. 
1993) showed reasonable agreement with the Langley 
nuclear, atomic and transport models. 

Experiments were performed at the 

Gessellschaft fur Schwerionenforschung (GSI) 

accelerator using beams of 12 C, 14 N, and 16 0 at energies 
of 674(±2)A MeV in which the transmitted flux of 
charge 5 to 8 were measured behind a water target of 
variable thickness (Schall, 1994; Schardt et al ., 1995; 
Sihver et al ., 1995). The measured transmitted fluence 
of all ions of the same charge as the beam (open circles) 
is shown in Fig. 16 along with the solution for the 
primary beam fluence (dashed curve) and the calculated 
fluence of all ions of charge equal to the initial beam 
(solid curve). It appears that the total absorption and 
neutron removal cross sections of the NUCFRG2 are 
reasonably correct (see also Golovchenko et al. 2002). 
The measured fluence of all ions with a single charge 
removed (open circles) is shown in Fig. 17. The 
NUCFRG2 code including alpha knockout processes 
shown as the first collision term (dashed) and the 
complete solution (solid) is in reasonable agreement 
with the experiments. The measured fluence for 
removal of two charge units from the initial beams of 
14 N and 16 0 (open circles) is shown in Fig. 18 with the 
NUCFRG2 results. The effects of alpha clustering is 
most apparent in the alpha knockout process for 16 0 
collisions (Figure 18, right panel). The carbon isotope 
distribution in highly nonelastic collisions are equally 
distributed between 12 C and 13 C as can be seen in 40 Ar 
fragmentation (Tull, 1990). In distinction, the 
fragmentation of 16 0 shows the single alpha knockout 
cross section to cause an excess of 12 C fragments being 
produced (Olson et al ., 1983, Cucinotta and Dubey 
1994, Wilson et al. 1998). The addition of the alpha 
knockout cross section leaving the 12 C core in the 
ground state to the NUCFRG2 nonelastic cross section 
(solid curve) brings good agreement with the GSI 
oxygen beam data. The carbon fragments produced on 
the first collision is discussed by Wilson et al. (1998). 
The only triple charge removal data is for the 16 0 ion 
beam as shown in Fig. 19 for which the NUCFRG2 is in 
reasonable agreement. 

High resolution electronic detectors have the 
added advantage of allowing the detailed testing of the 
nuclear and atomic interaction database as well as the 
computational models. The detector arrangement is 
shown in Fig. 20. The detectors upstream from the 
target are used to identify the incoming ion to insure it is 
in fact a beam particle and not a fragment produced 
along the beam pipe. Following the target are additional 
analyzing detectors to characterize the particles leaving 
the target. The detector responses are determined by the 
solution of equation (1) and are closely related to the 
straggling solutions in equation (8). Although multiple 
scattering was prominent in the older setup, such effects 



Target 


Fig. 20. Detector configuration (not to scale). All detectors are 
silicon; the trigger detector TR detector effectively defines the 
size of the usable beam spot. The respective thicknesses and 
approximate active areas are: TR - 300 pm x 300 mm 2 ; PSD1 
and PSD2 - 1mm x 1500 mm 2 ; d3mmU and d3mml-4 - 3mm x 
450 mm 2 . 



Fig. 21. Summed energy loss (MeV) in the detector pair d3mm3 
and d3mm4 for fragments produced by 1037 A MeV 56 Fe ions in 
10 g/cm 2 graphite-epoxy (50.92/49.08%). 



Fig. 22. Summed energy loss (MeV) in the detector pair d3mm3 
and d3mm4 for fragments produced by 1037 A MeV 56 Fe ions in 
5 g/cm 2 graphite-epoxy (50.92/49.08%). 



Fig. 23. Summed energy loss (MeV) in the detector pair d3mm3 
and d3mm4 for fragments produced by 1037 A MeV 56 Fe ions in 
2.54 cm Al. 
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are minimized in the present apparatus. In modeling the experiment described above, it was assumed that the 56 Fe 
beam extracted from the AGS at 1037 A MeV was nearly monoenergetic (the small inherent energy width had no 
influence on subsequent comparisons as the main width comes from the detectors and the target). The Green’s 
function solution (Tweed et al. 2005) was used to estimate the flux of the fragments entering the detector pair 
d3mm3 and d3mm4, which were treated as a single 6mm detector. The energy deposited by each fragment was then 
computed by means of equations (8) and used to evaluate the energy loss spectrum using the NUCFRG2 nuclear 
database (Wilson et al 1987b). The results (Tweed et al. 2003, Walker et al. 2005) are exhibited in Figs. 21 through 
23 where the computed energy loss spectrum (red curve) is compared with the experimental measurements (black 
curve). 

The right hand peak is from the surviving ion beam and fragments produced by neutron removal. The area 
associated with this peak is determined by the nuclear absorptive and neutron removal cross sections. The peak 
location and width are determined mainly by the range-energy relations and straggling process respectively. Clearly 
the comparison of the right most peak to the model is a strong test on these limited processes. The next peak to the 
left is for Mn fragments followed by Cr fragments and so on. These successive peaks are mainly governed by the 
fragmentation process and the area relates to the fragmentation cross section while the location and width is mainly 
determined by the energy down shift and energy dispersion on fragmentation. Hence the most accurate evaluation 
of the computational model lies in such comparisons. It appears that the Mn fragments are under-represented by the 
NUCFRG2 model for graphite/epoxy targets, which is consistent with earlier cross section measurements (Zeitlin et 
al. 1997). Overall, the fragmentation cross sections in Al seem better represented by NUCFRG2 than those for the 
graphite-epoxy composite. More recent results from the Quantum Multiple Scattering Fragmentation (QMSFRG) 
model may give some improvements (Cucinotta et al. 2003) over results in Figs. 21-23. An additional discrepancy 
is the energy downshift parameter for large mass removal especially for vanadium and lighter fragments. A further 
limitation of the model is for the lightest fragments (especially hydrogen and helium) that are produced with a 
broad spectrum resulting in a broad range of energy losses on the left margin of the measurements. This broad 
spectrum will be treated in later model development. The next emphasis will be on effects of multiple scattering 
and off-axis production components of ion fragmentation. Such effects will give improvements in the valley 
regions of the modeled results shown in Figs. 21 through 23. 

In order to further validate the model, it is desirable to compare with more than just the energy deposited in 
the detectors but to compare integrated quantities used in radiobiology to express the character of the beam. It is 
advantageous to make comparisons of the computational model with data that is derived from the fluences 
emerging from the target. The easiest and perhaps most important of these integrated quantities to compare with is 
the track average LET, (LET) trk . The (LET) trk can be computed from the formula 

(LET) trk = I, J dQ dE Lj(E) <j)j(x,Q, E)/I, j dQ dE <()j(x,Q,E) (37) 

where <bj(x,Q,E) is the flux and LEE) is the LET of ~ , „ . . , , ,, . . , , , „ 

TJ J Table 1: Track Averaged LET from Model and Experiments. 

they particle type at energy E. In order to be able to Angular factor for z=l (2) is 0.02 (0.08) 


make comparisons with experimental data, we will 


Depth 

<LET> trk 

<LET> trk 

take S/E) to be that in water for particle j at energy 

Material 

(g/cm A 2) 

model* 

experiments 

E\ i.e., the stopping power for water. 

Carbon 

3.9 

125.3 

127.0 

The light ions, Zj = 1,2 have a strong 

Aluminum 

7 

127.3 

125.4 

angular dependence and, as a consequence of the 

Lead 

3.6 

148.2 

145.8 

straightahead approximation currently used in the 

Polyethylene 

10 

91.3 

91.4 

model, there is a marked over prediction in their 

Graphite-Epoxy 

5 

116.3 

121.3 

fluences. If uncorrected, the predictions for the 

Graphite-Epoxy 

10 

94.8 

98.5 


(LET) trk would be drastically lower than 

experiments, due to the artificial inflation of the denominator in equation (37). It is possible to correct for this by 
scaling Z 7 = 1,2 by an appropriate weight factor for each of the light ions, thereby reducing their number and 
dampening the effect. The weight factors are 0.02 for Zj = 1 and 0.08 for Z 7 = 2. These corrections are used for all 
targets and 56 Fe beams (Walker et al. 2005). Our ability to improve on these estimates will be limited by the lack of 
a supporting nuclear physics program with improvements as time available basis. 

Table 1 shows (LET) trk computed from the model with the corrected fluences and preliminary results from 
experiments done at Brookhaven National Laboratory Alternating Gradient Synchrotron (BNL AGS) for various 
targets and thickness. All of the experiments were carried out with a 1 A GeV 56 Fe beam. In all cases, there is 
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reasonable agreement between the predictions of the model (± 5 percent) and the experimental results with the 
greatest error occurring in the results for graphite-epoxy. 

Despite the clear advantage of direct comparison to detector output and its relation to atomic and nuclear 
processes, it has been the decision of the NASA Space Radiation Shielding Program to work with reduced 
quantities such as those given in Table 1. Recently, a benchmark was defined by 1 A GeV 56 Fe ion beams in 
polyethylene, polymethyl metachrylate (PMMA), aluminum, and lead at several depths. The benchmark quantities 
requested are for shielding materials program are: 

Dose (unit 56 Fe fluence): 

D(x) = Klj / dG dE L/E) (/)j(x,G,E) 


Track averaged LET: 

(LET) trk = Zj JdG dE L/E) (j)j(x,G,E) 
/ZjfdGdE c/)j(x,G,E) 


Dose averaged LET: 

(LET) dose = Ij Ida dE L/E) 2 <j)j(x, Q, E) 
/Ij Ida dE L/E ) <pj(x,a,E) 


Fraction of Fe-ions: 


Fsurviving (/) JdG dE (j)j(x, G,E) 


/ incident ion fluence (unity) 


For the targets of polyethylene, poly methyl 
metachrylate (PMMA), aluminum, and lead 
over the thickness range of 5 to 30 g/cm 2 . In 
addition to the quantities requested, fragment 
spectra in PMMA as 

Fragment spectrum at 23 g/cm 2 PMMA: 

Fj(x,E) = JdG (j)j(x,G,E) 


/ incident ion fluence (unity) 

where K is conversion from MeV/g to nGy 
(K=1.602). In that the full angular dependence 
is lacking in the present computational model 
(Tweed et al. 2005), we have fit a 
renormalization constant to the light ion data of 
Zeitlin, Miller, and Heilbronn with resulting 
values of 0.02 (Z=l) and 0.08 (Z=2) as 
discussed by Walker et al. (2005). The dose, 
track and dose averaged LET, and fraction of surviving Fe-ions are given in Table 2. The iron ion range, R 0 , for 
each target material is given in the table. The fragment spectra penetrating 23 g/cm 2 of PMMA are given in Fig. 
24. The energy spectra of the penetrating iron ions and iron fragments are shown in greater detail in Fig. 25. The 
shift of the iron fragment spectra to lower energies relative to the surviving beam particles resulting from loss of 
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Table 2. Benchmark parameters for 1 A GeV iron ions in the four 
requested materials. 

Polyethylene (Rq= 25.3 g/cm 2 ) 

Depth 

(g/cm A 2) 

<LET> trk 

(keV/pm) 

<LET> dose 

(keV/|im) 

Dose 

(nGy) 

Fraction Surviving 
Fe 56 Fe 55 ' 54 Total 

5 

112.9 

142.4 

180.8 

0.50 

0.030 

0.53 

10 

92.3 

139.2 

147.9 

0.25 

0.032 

0.29 

15 

82.3 

143.4 

131.9 

0.13 

0.025 

0.15 

20 

81.1 

163.8 

129.9 

0.067 

0.019 

0.085 

25 

98.5 

274.9 

157.8 

0.035 

0.0076 

0.043 

30 

33.0 

295.7 

52.8 

0 

0 

0 

PMMA (Ro= 26.1 g/cm 2 ) 

Depth 

(g/cm A 2) 

<LET> trk 

(keV/pm) 

<LET> dose 

(keV/|im) 

Dose 

(nGy) 

Fraction Surviving 
Fe 56 Fe 5554 Total 

5 

116.1 

143.7 

186.1 

0.54 

0.029 

0.57 

10 

96.3 

141.0 

154.1 

0.29 

0.032 

0.32 

15 

86.1 

144.8 

138.0 

0.16 

0.027 

0.19 

20 

84.7 

162.5 

135.7 

0.088 

0.021 

0.11 

25 

116.9 

301.7 

187.2 

0.050 

0.015 

0.065 

30 

42.3 

334.3 

67.8 

0 

0 

0 

Aluminum (R 0 = 33.9 g/cm 2 ) 

Depth 

(g/cm A 2) 

<LET> trk 

(keV/pm) 

<LET> dose 

(keV/|im) 

Dose 

(nGy) 

Fraction Surviving 
Fe 56 Fe 5554 Total 

5 

133.3 

150.1 

213.5 

0.78 

0.0093 

0.79 

10 

119.2 

151.3 

191.0 

0.61 

0.015 

0.63 

15 

109.4 

154.6 

175.2 

0.48 

0.017 

0.50 

20 

104.5 

162.3 

166.7 

0.38 

0.018 

0.40 

25 

104.5 

179.1 

167.4 

0.30 

0.018 

0.32 

30 

120.3 

229.9 

192.7 

0.24 

0.017 

0.25 

Lead (Rq= 50.5 g/cm 2 ) 

Depth 

(g/cm A 2) 

<LET> trk 

(keV/pm) 

<LET> dose 

(keV/|im) 

Dose 

(nGy) 

Fraction Surviving 
Fe 56 Fe 5554 Total 

5 

147.2 

152.0 

235.8 

0.93 

0.0015 

0.93 

10 

144.0 

154.3 

230.6 

0.86 

0.0028 

0.87 

15 

141.5 

157.3 

226.7 

0.80 

0.0040 

0.81 

20 

139.9 

161.5 

224.2 

0.75 

0.0050 

0.75 

25 

139.9 

167.3 

224.1 

0.69 

0.0058 

0.70 

30 

141.9 

175.8 

227.3 

0.65 

0.0065 

0.65 




inertia is clearly seen in the graph. Similar observations were made in the neon beam experiments of Schimmerling 
et al. (Shavers et al. 1993). 

SPACEFLIGHT VALIDATION 
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The present benchmark is a useful exercise for 
comparison of computational procedures and 
atomic/nuclear database. The first year of the 
current project has been mainly focused on 
advancing computational procedures (Tweed et al. 

2005), validation of those procedures (Walker et 
al. 2005), and providing a sensitivity analysis on 
the 14,365 nuclear parameters required for space 
radiation simulations in which 170 isotopes are 
represented (Heinbockel et al. 2005). This has left 
little time for improving the physical description 
of the transport process that is further hampered 
by the need of a supporting basic nuclear physics 
research. The main progress is towards 
implementing multiple Coulomb scattering into 
the formalism to be followed by off-axis nuclear 
fragmentation components with a simplified 
nuclear model. The first issue to be addressed 
prior to flight validation is the sensitivity of the 
transport solutions on the nuclear parameters. The 
theory of nuclear absorption cross sections is well 
known and new experimental assessments 
consistently show these cross sections to be within 
5 percent or better as new measurements are made 
(Tripathi et al. 1998,1999, Golovchenko et al. 

2002). However, the distribution of the reaction 
over various reaction channels are less accurately 
known but experimental knowledge is improving. 

Since the transport of ions involves the integration 
of well-understood atomic processes and less 
well-understood nuclear processes, a parametric 
sensitivity study to understand the role of nuclear 
cross section uncertainties is in order. In such a 
sensitivity analysis, the nuclear absorption cross 
sections are assumed to be well understood and 
left unchanged while the specific channel cross 
section sensitivity will be ascertained. 

Although the biological response to HZE 
exposure is very uncertain, dose equivalent is still 
assumed a first order approximant of biological 
injury (NCRP 2001) and will be used herein to 

study fragmentation cross section sensitivity. The dose equivalent is given as 
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Fig. 24. Iron ion fragment spectra behind 23 g/cm 2 of PMMA. 
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Fig. 25. Iron ion spectra after 23 g/cm 2 of PMMA. 


H(x) = KIj JdQ dE L/E) Q[ L/E) J (fj(x, Q, E) 


(38) 


where Q[Lj(E) ] is the LET dependent quality factor. The particle fields (f)j(x, Q, E) contain uncollided particles and 
secondary particles which depend on specific reaction cross sections. Of particular importance are the HZE 
particles with their large LET and associated quality factors providing large contributions to dose equivalent. For 
example, the Fe ions of the galactic cosmic rays are known to contribute about one third of the total dose equivalent 
inside space structures (Wilson et al. 1991). We explicitly display the dependence of dose equivalent on the 
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fragmentation cross sections while suppressing the dependence on x and expand to first order about nominal cross 
sections given by NUCFRG2 as follows 

H(a jk + £ jk ) = H(o ik ) + d^HtOpJ e jk (39) 

where £ Jk represents the uncertainty in the NUCFRG2 fragmentation cross section Oj k . Note, the term () ol iJ I ( <J !k ) 
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represents the sensitivity of H to the variations in the cross section parameter o jk . Equation (39) has been used as a 
basis for a study of the uncertainty propagation into the solution as a function of shielding. 

The sensitivity coefficient d^k^Ojk) is a complicated function of the 14,365 cross sections Oj k and the 
particles present at the boundary. These are evaluated numerically by finite differencing. It has long been known 
(Wilson et al. 1997) that only seven ions give over 90 percent of the contributions to dose equivalent: H, He, C, O, 
Mg, Si, and Fe. Clearly, these ions and their secondaries will play a major role. The sensitivity coefficients 
dq k H(Oj k )XGj k are shown in Fig. 26. The results shown are for aluminum shielding but other shield compositions are 
qualitatively similar. The contributions from a single ion type at the boundary is indicated by the sum of 
contributions of all fragments from the same type of incident ion at the boundary. These are also shown in Fig. 26 
by separate panel as the cumulative error for each incident ion type. The essential role of the seven important ions 
before identified is clear. In the case of iron, the fragments for charge removal up to about 10 are seen as the most 
important contributors to dose equivalent even up to depths of 40 g/cm 2 . This is due to their abundance and large 
charge and quality factors. At the larger penetration depths, the light ion fragments from incident iron ions play an 
increasing role but are dominated by the heavy fragments even at the large depths shown in the figure. The ions of 
Mg and Si still have significant contributions from fragments with few to several charges removed but the light ions 
are relatively more important in this case than in the case of iron ions. The ions below Mg have relatively large 
contributions from the light ion fragments they produce especially at the large depths. It is fortunate that only 
relatively few fragments of each of relative few incident ions need to be known with a degree of accuracy. We now 
examine what impact uncertainty in fragmentation cross sections and what level of certainty is required for an 
adequate solution of the shielding problem. 

A cross section uncertainty effects study on the propagated errors in polyethylene, aluminum and copper 
were performed using the sensitivity data like that in Fig. 26. The propagated error associated with each 
fragmentation cross section uncertainty is p xd q j k H(G jk )XG jk where p is a random sequence representing the 
distribution of propagated uncertainty and chosen from a normal distribution of zero mean and 0.75 standard 
deviation (a very conservative assumption). The results on the propagated error in the three materials are shown in 
Fig. 27 as 1 -a error bars for the propagated uncertainties (Heinbockel et al. 2005). The percent rms propagated 
error is given in Table 3. 


Table 3. Propagated percent rms error due to 75 percent 
uncertainty in ion fragmentation cross sections. 


Depth, g/cm 2 

Polyethylene 

Aluminum 

Copper 

10 

3.49 

1.18 

0.70 

20 

3.58 

1.85 

0.96 

30 

3.18 

2.02 

1.39 

40 

2.00 

1.86 

1.53 

Also shown 

in Fig. 27 are 

solutions with all of 



Ojk set to zero for comparison. Note, the zero bias 
and non-systematic assumption leave little 
uncertainty in the propagated error. We have done 
systematic error studies which expand the error bars 
in an unsymmetrical way but still the width of the 
error bars are minor. It is clear from this result that 
any added gains in accuracy in fragmentation 
parameters beyond the NUCFRG2 model will 
provide minimal improvements to space radiation 
shield evaluation. However, the continued 
development of transport procedures and 
atomic/nuclear database has value in beam 
characterization for the radiobiology program. The 
above error propagation analysis will help to make 
the results of validation with spaceflight data better 
understood. 

There are two principal flight platforms used in the spaceflight validation of the shield analysis codes: the 
space shuttle and the International Space Station. There is also flight data from space station Mir but the 
knowledge of the material arrangement on Mir is limited and therefore not used in quantitative analysis. The 
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Fig. 27. Nominal dose equivalent and cross section 
uncertainty comparison with production cross-sections equal 
to zero (Ts = 0). 


instrumentation has varied over the time of 
measurement. Although some preliminary 
data has been analyzed from the Odyssey 
Mars Radiation Environment Experiment 
(MARIE) instrument in Mars orbit, most 
validation data has been gathered in low Earth 
orbit (LEO). The space shuttle has well 
defined geometry models consisting mostly of 
aluminum alloy 2219 although loading on 
specific missions which also has well 
characterized geometry is important in 
specific mission data (for example, STS-61A 
the first German Spacelab mission Dl). An 
adequate ISS shield model is developing with 
the validation of the ISS 7 A configuration 
(Hugger et al. 2003a, b). A historic problem 
with LEO validation is the uncertainty in the 
trapped particles (mainly protons). Available 
proton models denoted AP8 MIN and AP8 
MAX were developed to represent the 
omnidirectional proton flux near solar 
minimum (circa 1965) and solar maximum 
(circa 1970). This model uncertainty was 
partly corrected in STS and ISS dosimetry 
analysis by adjustments to the trapped flux 
according to thermoluminescence detector 
(TLD) readouts. The degree to which AP8 
uncertainty enters depends on the orbit 
inclination and altitude. 

The two primary limitations in the 
LEO trapped environmental models AP8 
MIN&MAX as discussed by Wilson et al. (2002d) is the assumption that the trapped particles are isotropic 
(resulting from the omni-directional fluence description) and poor representation of the dynamic behavior. These 
omnidirectional models have been relatively successful in describing the radiation environment aboard the highly 
maneuverable Shuttle wherein anisotropies tend to be averaged. Such models will not be adequate in the formation 
flying of the ISS mainly oriented in the local horizontal plane along the velocity vector. We have developed a 
dynamic/anisotropic trapped proton environmental for future use in LEO shield design and operations. These 
environmental models are placed in a suitable form for evaluation of the incident radiation on the bounding surface 
of the six-degree of freedom motion of an orbiting spacecraft for shield evaluation (Wilson et al. 2005b). To test 
the dynamic behavior, we use shuttle TLD data available from 1981 to the present giving good coverage for over 
two solar cycles. A sample of shuttle TLD measurements and the computational model is given in Table 4. A 
compilation by Badhwar (Chapter 21 of Wilson et al. 1997) of calculated and observed dose rates in shuttle flight is 
given in Fig. 28. According to Badhwar (2002), the rms error of both observed and calculated dose rate and dose 
equivalent rate is less than 15 percent. 

Many detector systems are designed to capture some physical quantity such as LET but do so only within 
the limitations of the detector. Generally we have approached this problem by not comparing to the quantity that 
the detector is meant to measure but rather to evaluate the detector response to its environment as estimated using 
shielding models of the platform in which the measurement was made. One such example is the nuclear track 
detector (CR-39) measurements on the Dl mission of STS-61 A. The results of computational models with detector 
response functions which depend on details of post-flight processing of the exposed detectors as estimated by Shinn 
et al. (1995) is shown in Fig. 29. Many of the target fragments at the highest LET were not seen by the CR-39 
plastic as they were etched away in the track development process but the predicted CR-39 and measured CR-39 
responses are good agreement. 

A great boon to model validation was the development of the particle detector telescope spectrometer for 
use on the shuttle (Badhwar et al. 1992, 1994). This high-resolution detector allows detailed testing of spectral 
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Table 4. Comparison of present model with Shuttle flight data. 


Flight 

Date 

DRNM* 

DLOC 

TLDf 

(|xGy/d) 

Calc. 

(|xGy/d) 

STS-41A 

11/83 

6421 

3 

64.6 

60.8 

STS-31 

4/90 

5701 

1 

2141 

2093 

STS-62 

3/94 

6771 

1 

94.3 

90.7 

STS-65 

7/94 

6822 

2 

28.3 

25.7 

STS-82 

2/97 

7074 

1 

2978 

3014 

STS-92 

10/00 

6417 

2 

165.9 

172.9 


* Deep River Neutron Monitor count rate f GCR corrected TLD -100 data 



Fig. 28. Observed dose rate versus calculated dose rate (Badhwar, 
Chapter 21 of Wilson et al. 1997). 




components transmitted through or produced in the shuttle hull. Only the more intense components can be detected 
with sufficient statistical accuracy on the relative short shuttle flights. Figs. 30 and 31 show particles produced in 
the shuttle hull in the high latitude regions of that particular flight (STS-48) where galactic cosmic rays are present 
in sufficient intensity although the geomagnetic cutoffs are still relatively high. The particles in the figures 
represent the secondaries produced in the shuttle hull by much higher energy ions not present in the measurements 
or the graphs. The older version of HZETRN did not include the direct cluster knockout particles and the figures 
show mainly decay particles of target nuclei of the shuttle hull material struck by the incoming ions. The revised 
and current HZETRN include these knockout processes and show excellent agreement with measurements 



Fig. 29. CR-39 detector response (triangles) compared to 
predicted CR-39 response and LET spectra during the 
German Spacelab mission D1 on STS-61A. 



Fig. 30. Comparison of calculated and measured proton 
spectra on STS-48. 
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Fig. 32. Galactic cosmic ray induced lineal energy 
spectrum during STS-57. 
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Fig. 33. Trapped proton induced lineal energy 
spectrum during STS-57. 



Fig. 3 1 . Comparison of calculated and measured deuteron 
spectra on STS-48. 


Fig. 34. Comparison of HZETRN with multigroup 
methods of Clowdsley with the old and new JSC 
analysis of STS-36 measurements. 
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and the current HZETRN include these knockout processes 
showing excellent agreement with the measurements and 
excellent validation of the light ion production processes. 

The Tissue Equivalent Proportional Counter (TEPC) 
provides spectral distributions of energy deposited within a 
small volume of gas to approximate the quality of the 
radiation absorbed (usually represented by LET). The TEPC 
spectrum is the defined “lineal energy spectrum” as it does 
not exactly but does somewhat approximates an LET 
distribution (see Badhwar et al. 1994). The usefulness of 
this instrument had been limited for many years for lack of a 
useful model of its response to the space radiations. Shinn 
et al. (2001) have developed a model for the response of the 
TEPC used by Badhwar for an estimate of dose equivalent 
in spaceflight. By developing a model for the response of 
the TEPC, Shinn et al. hoped to make this instrument more 
useful in testing shielding models using spaceflight data 
since flying the TEPC became standard practice on shuttle 
flights for the last ten years. Further usefulness of the TEPC 
results from time serial analysis of the measured data to 
separate the trapped environment from galactic cosmic ray 
contributions allowing a separate test on shielding against 
trapped protons and shielding against galactic cosmic rays. 
Analysis of the environment induced in shuttle flight 
STS-57 by galactic cosmic rays and by trapped protons 
as measured and calculated using HZETRN are shown 
in Figs. 32 and 33. It is expected that the TEPC 
estimates of dose equivalent will be in good agreement 
with the computational models as Badhwar has 
determined. 

Limited measurements with a Bonner sphere 
setup using foil activation methods have been made 
(Keith et al. 1992) with re-analysis of Keith’s data by 
Badhwar et al. (2001). The old and new analysis of the 
Bonner sphere data are shown in Fig. 34 with the 
revised HZETRN with multigroup methods of 
Clowdsley et al. (2000). Clearly the computational 
model is within the uncertainty of the data analysis. 
Although the results in Fig. 34 are useful in evaluation 
of computational models, an improved approach would 
be to calculate the response of the four Bonner spheres 
from computational models with direct comparison with 
the foil count rates of the measured sphere data. This 
would eliminate the need for assumptions as to the 

Table 5. Calculated and measured Quality Factors 


Organ 

Qcalc 

Qmeas 

Brain 

2.10 

1.8±0.1 

Colon 

2.17 

2.0±0.2 

Heart 

2.15 

2.0 ±0.4 

Stomach 

2.2 

1.7 ±0.3 

Thyroid 

1.89 

1.7 ±0.2 

Skin breast 

2.10 

1.8 ±0.2 

Skin abc iomen2. 1 1 

1.9 ±0.2 



Fig. 35. The JSC phantom torso as flown on STS-91. 
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Fig. 36. Apparatus for shield attenuation studies during 
shuttle flight. 


STS Data in Payload Bay 



Fig. 37. Experimentally measured and calculated attenuation 
curves in aluminum and polyethylene in shuttle flight. 
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Fig. 38. The ISS 2A.2B Configuration CAD Model. 


neutron spectral components as input to the Bonner sphere analysis methods giving rise to the uncertainties seen in 
the figure. 

Applying radiation shield constraints to the space 
radiation shield design process requires not only the 
estimation of the internal spacecraft environment but how 
that environment is transported in human tissues of a 
characteristic human shape. To examine this problem a 
human skeleton (head and torso) is covered with a tissue 
equivalent plastic and sectioned to allow instrumentation for 
in-flight measurements (see Fig. 35). The torso phantom 
was CT scanned to define the geometry for use as a 
computational model. The instrumentation used TLD, 
nuclear track foils and TEPCs to allow the mapping of the 
dose and dose equivalent within the phantom torso model. 

The instrumented torso phantom was flown on STS-91 
followed by data analysis and HZETRN validation. The 
results of the measurements expressed herein as the average 
quality factor is given in Table 5. For the listed organs, the 
calculated dose was within 12 percent of the measurements 
except for the stomach that was within 19 percent. The 
cause of the differences was in large part due to need to 
correct the galactic cosmic ray under-response of the TLDs 
(Cucinotta personal communication, 2004) and the corrected 
data was all within 12 percent of the computational models. 

The usefulness of the TEPC in evaluating dose 
equivalent in the space environment being established, 

Badhwar provided a test rig of five TEPCs (Fig. 36) to study 
dose equivalent attenuation in the LEO environment. The 

bare detector allows characterization of the ambient field within the shuttle while the shielded detectors sample the 
modified environment due the shielding material. A study of LEO radiation attenuation in polyethylene was 
conducted in the shuttle bay with the doors open on STS-81 and are compared to a similar study of aluminum 
shielding on STS-89. Note that a direct comparison of aluminum shielding effectiveness to that of polyethylene 
cannot be made as the two flights are in radically different environments. STS-81 is in a high inclination orbit 
(51.6°) with large contributions from galactic cosmic radiation, and STS-89 is in a low inclination orbit (28.5°) 
where the environment is dominated by trapped protons. The results are compared with HZETRN model 
calculations in Fig. 37. One TEPC failed on STS-89 but it is still clear that good agreement was found between the 
HZETRN modeled and the measured attenuation curves. This is also true for the STS-81 results for polyethylene 
shields as shown in the figure. It is clear that polyethylene has superior attenuation characteristics even in 
environments with large contributions from galactic cosmic radiation (STS-81). As a result, polyethylene is the 
material of choice for augmenting the ISS shield. Most important is the ability of HZETRN to evaluate correctly 
the dose equivalent from all components important to human rated systems and provide estimates of ISS 



Fig. 39. The ISS 7 A Configuration CAD Model. 


Table 6. TLD dosimeter locations. 


Location 

Description 

Labi 

Starboard side standoff between starboard and ceiling 

Lab2 

Near the Lab Window 

Lab3 

Port side standoff between Deck and Port rack 

Lab4 

Forward closeout panel starboard upper corner 

Lab5 

Sleep Station 

SMI 

Inside Port SM Crew Quarters, Outboard wall aft upper corner 

SM2 

Inside Starboard SM Crew Quarters - Outboard wall aft upper corner 

SM3 

Panel 339, aft section behind treadmill- Upper center part of the panel 

SM4 

Panel 32 7, overhead, forward of treadmill 

SM5 

Panel 307 overhead, near center 

SM6 

Window #14, Transfer Compartment Adapter Section, Starboard Nadir Quadrant 
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augmentation requirements and value of the benefits. 

The omnidirectional models of the trapped 
environment have been useful for highly maneuverable 
spacecraft such as the shuttle. For spacecraft like the ISS 
with its relatively fixed flight attitude in the local vertical 
reference frame so there is no angle averaging process, the 
direction of arrival of radiation on the spacecraft exterior 
is expected to play an important role. The ISS analysis 
started with use of the dynamic omnidirectional physics 
based models connected to the conventional AP8 
MIN&MAX database environments for 1965 and 1970 
respectively (Wilson et al. 2002d) 

To create the shielding model, the ISS was re- 
created using the computer aided design package I-DEAS. 
The ISS is an on-going endeavor, which has taken the ISS 
through many stages. To accommodate the changing 
space station, the CAD shielding model representing the 
ISS must be re-configurable as well. To do this, the 
individual components were modeled separately and then 
combined before analysis into the necessary configuration. 
All of the components from the first ISS configuration to 
the current ISS 11A configuration have been modeled. 
The Components are the Service Module, the Lab Module, 
the Node 1, the FGB, the 3 Pressurized Mating Adaptors, 
the Airlock, the Docking Compartment, the Soyuz, and the 
Truss. For the purpose of this paper, the ISS shielding 
model needed to represent all of the ISS configurations 
from the 2A.2B through the 7A. The 2A.2B and 7A 
configurations are shown in Figs. 38 and 39. 

Each ISS component is comprised of multiple 
parts. In modeling the parts of each component, the focus 
was on simplicity while trying to maintain the accuracy. 
Hence, the part’s solid angle and thickness relative to the 
target point had to be preserved (Qualls et al., 1997). To 
do this, first it was necessary to find out the masses and 
dimensions of all the objects in an ISS component. For a 
better description of how individual parts of ISS 
components were modeled for the ISS shielding model, 
see Hugger, (2003a). 

The initial dose calculation was performed only 
for the omni-directional trapped radiation case (Wilson et 
al. 2002d) for the dosimetry locations in Table 6. The 
anisotropy of the trapped radiation environment was later 
introduced using the angular dependence in the South 
Atlantic Anomaly averaged over the ascending and 
descending node crossings. The GCR environment 
includes the Earth’s shadow in both simulation modes. 
The six dosimeter locations in the Service Module are 
shown in Fig. 40. Five of the dose locations (SM-1 to 
SM-5) are well shielded in all directions. The location 
SM-6 however is in the pressurized docking adapter with a 
significant fraction of the modeled shield below 10 g/cm 2 
(Hugger et al. 2003b). The dose along each ray was then 
calculated by interpolating the amount of shielding within 
dose versus depth data with specified amount of backing. 



Fig. 40. Dosimeter locations in the ISS Service Module. 


Table 7. Time periods and configurations. 

Period 

Dates 

Altitude (km) 

Time 1 

09/08/00-12/11/00 

385.75 

Time 2 

12/01/00-05/01/01 

381.04 

Time 3 

04/19/01-07/24/01 

390.54 

Time 4 

07/12/01-08/22/01 

397.02 

2A2B 

09/08/01-10/11/01 

373.2 

3A 

10/11/00-10/31/00 

389.9 

2R 

10/31/00-11/30/00 

402.8 

4A 

11/30/00-02/07/01 

393.5 

5A 

02/07/01-04/19/01 

398.9 

6A 

04/19/01-07/12/01 

401.8 

7A 

07/12/01-08/22/01 

389.6 

7A1 

08/10/01-08/22/01 

396.8 


Table 8. Days for each configuration in a Time Period. 


Configuration 

Period 1 

Period 2 

Period 3 

Period 4 

2A2B 

34 

- 

- 

- 

3A 

20 

- 

- 

- 

2R 

29 

- 

- 

- 

4A 

10 

69 

- 

- 

5A 

- 

71 

- 

- 

6A 

- 

11 

84 

- 

7A 

- 

- 

12 

29 

7A1 

- 

- 

- 

11 


SM Dose Comparison 



Fig. 41. Comparison of Methods 1 & 2 with JSC flight 
data measured in the service module. 
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Anisotropic vs Isotropic Environment 



Anisotropic 


Isotropic 


□ Method 1 ■ Method 2 □ Measured 


Fig. 42. Anisotropic effects on SM-6 dosimetry. 


This process was repeated for all the rays in the thickness file. The total dose was calculated by summing all of the 
directional doses and then dividing by the total number of directions (depends on normalization). 

The measurements took place over four separate time periods using thermo-luminescence detectors (TLD- 
100) located throughout the ISS. The measurements were performed over the course of nearly 11 months. The 
dates and time periods for the different measurements are listed in Table 7 with the associated ISS configurations. 
During these measurements, the ISS was constantly changing configuration, altitudes, and the total number of 
dosimeters used was changing as well. For the first testing period, there were four TLD’s measuring data in the 
Service Module. The second testing period had 6 TLD’s operating in the SM. The third and fourth testing periods 
utilized all 1 1 TLD’s in the Service and the Laboratory modules. The TLD locations are given in Table 6. 

To accurately model the environment for each set of dosimeter data, the altitude and time period were 
required. The periods labeled Time 1 through Time 4 are periods for accumulating dosimeter data. The altitudes 
listed in this case are the averages over the as flown five-minute trajectory data. The ISS configurations with 
associated dates and altitudes were provided by the Johnson Space Center. For each configuration of the ISS, an 
average altitude was calculated and an 
environment was generated. There were a total 
of 8 environments for eight ISS configurations 
created, also listed in Table 8. The number of 
days of operation for each configuration within 
the listed time periods are given in Table 7. For 
this paper, we will concentrate on time period 4 
for which the configuration nearest the current 
11A and use various methods for comparing 
with the TLD dosimeter data. No correction was 
made for under-response of the TLD to highly 
ionizing particles, resulting in -10% 
underestimate of absorbed dose in water. 

The first evaluation method (Method 1) 
uses the average altitude for the time period 4 
but allows evaluation under the two 
configurations 7A and 7A1. A second method 
(Method 2) used the altitude apogee and perigee 
data provided by Johnson Space Center applied 
to each configuration as given in Table 8 for 7A 
and 7A1. In each case, the dose was evaluated 
at each of the six dosimeter points in the Service 
Module (SM-1 to SM-6) and for five locations 
in the Laboratory Module (Lab-1 to Lab-5). 

Absorbed dose in water calculated for 
the SM locations using Method 1 & 2 are shown 
in Fig. 41 in comparison with the dosimetry data 
measured in the same time period. The TLD 
dose is used without the small correction for 
under-response to high LET events. The only 
difference between methods 1 & 2 is the 
altitudes used in generating environmental data 
for the ISS configurations in Table 8. The 
differences can then be seen as our uncertainty 
in the actual flight altitudes, and it is highly 
desirable in further analysis to use as flown 
flight data. There is reasonable agreement 
between measurements (generally better than 20 
percent) and calculations except for the location 
SM-6 with its unusually thin shielding. Since 
the calculation was made using only omni- 
directional fluence data it appears that particle 
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Fig. 43. Comparison of Methods 1 & 2 with JSC flight data 
measured in the laboratory module. 



Fig. 44. Directional dose distribution predicted for ISS in passing 
through South Atlantic Anomaly. 
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arrival within the thin sections at SM-6 are overestimated by over 45 percent. It is well known that the largest 
fraction of space exposure comes from the least shielded region. 

To test the effects of anisotropy on the SM results, we implemented the directional dependent flux factor as 
given by Heckman and Nakano (1963). The directional intensity is then defined for direction of arrival defined by 
the local pitch angle 0 and azimuth X as follows 

j( Q,A) = J 4jt C exp[-( Q-n/2) 2 /2 Gq ] exp[r g cos(I) cos(X) // i s ] (40) 

where r g is the particle gyro-radius about the field line, I is the dip angle, h s is the atmospheric scale height, a e is the 
width of the pitch angle distribution, and C is a normalization constant as discussed by Kern (1994). The angular 
distribution is averaged over the South Atlantic Anomaly crossings for both ascending and descending nodes and 
the dose at SM-6 re-evaluated. The results are shown in Fig. 42 using both methods 1 and 2 for evaluation of the 
flux environment. It is clear that anisotropic effects are especially important for this location in the SM. 

A preliminary comparison with the Laboratory dosimetry has been made for the five locations Lab-1 to 
Lab-5 (see Table 6 for locations) during period 4. In this comparison, it is assumed that the trapped radiation is 
isotropic. The first attempt to compare with the Lab data left two storage racks empty due to lack of data at the 
time gives mixed results. Addition of 100 pounds of ill-defined “equipment” provides the results in Fig. 43. The 
agreement is reasonable for this first comparison (all of the data points are within 25%) but treating the anisotropy 
and the as flown trajectory data may lead to further improvement as it did in the case of the Service Module. 

The current status of LEO validation of shield modeling of ISS is encouraging. Clearly, the continued 
development of an understanding of ISS shield distribution using TLD measurements will require new software 
able to follow the six degrees of freedom trajectories in environmental evaluation using a fully dynamic and 
anisotropic trapped particle environment model (Wilson 
et al. 2005b). 

The effects of angular incidence of the LEO 
environment are shown in Fig. 44 on the directional 
dose sphere evaluated in the US Habitat module model. 

On the spheres in Fig. 44, one can clearly see the Earth 
shadow effects on the GCR arrival as the black region 
on the bottom of the sphere. The square figure on the 
front of the sphere results from the low shielding on the 
space end of the hallway. The diagonal slash is due to 
the arrival of protons in passage through the South 
Atlantic Anomaly (SAA). Clearly these features will 
have a strong influence on how the ISS shielding affects 
dosimetric measurements. One can not over emphasize 
the importance of these LEO studies as a test bed for 
shield design software, as test area for long duration 
human spaceflight, and as a staging area for future deep 
space missions. 

The current development has a two-fold 
purpose. First is the evaluation of the ISS shielding 
model using the ISS dosimeters. This first step is 
required to validate our understanding on the relation of 
the ISS distribution of materials on the internal radiation 
environment. Second is the evaluation of the best means 
of augmenting the Service Module Crew Quarters to 
lower the exposure of the crew during ISS operations. 

In both cases, we anticipate that the angular dependence 
of the arriving radiation is an important factor. 

The importance of angular dependence on ISS 
dosimeter results has already been demonstrated in some 
regions of the Service Module as shown for the TLD 
readings in the Transfer Compartment Adapter Section 
(SM06 dosimeter location, see Figs. 40 & 45) as shown 



Fig. 45. Cutaway view of Service Module showing 
dosimeter locations. 



Fig, 46. Cutaway view of ISS 1 1 A. 
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Fig. 47. Cosmic ray transition curves above 300 A 
MeV measured by Webber and Ormes (1967). 


in Fig. 42. In the figure, we show two methods of evaluating the 
environment on the basis of conflicting altitude descriptions for 
both assumed anisotropic and isotropic arrival. Method 1 uses 
the average altitude over the measurement period and Method 2 
uses the apogee and perigee evaluated separately and averaged. 

Independent of the remaining trajectory issues, it is clear that the 
anisotropy description is fundamental to understanding the ISS 
dosimeter data. The current code development will allow not 
only the inclusion of as flown ISS 3D trajectory data but the 
inclusion of limit cycle and maneuvers in the environment 
generation to increase the fidelity of future ISS shield studies. 

This is critical to evaluation of the ISS 1 1 A shield model shown 
in Figs. 45 and 46, the future augmentation of the ISS shield, 
and evaluation of the demonstration flight of the CEV in 2008, 
and the first uncrewed CEV flight but fully instrumented flight 
in 2011. 

In addition to space flight measurements, the transistion 
of galactic and solar cosmic rays in the Earth’s atmosphere has 
been studied for many decades. Early limitations of state-of-the- 
art nuclear database generators was a result of comparisons with 
high altitude balloon flight data gained with gas Ceremkov 

counters (Webber and Ormes 1967). In these experiments the ion flux of lithium and beyond was measured in the 
range of 200 A MeV to 5 A GeV with results shown in Fig. 47 in comparison to HZETRN solutions with the 
original Bowman, Swiatecki and Tsang nuclear model (Bowman et al. 1973) and the modifications of the nuclear 
excitation energy model by Wilson et al. (1987b). These comparisons were essential to the development of the 
NUCFRG model. More recently Blattnig, Norbury, and 
Cucinotta have continued development of pion production 
models and recently Blattnig and Norbury (Blattnig et al. 2004) 
have developed a version of HZETRN with meson production 
and transport and made comparison with atmospheric air shower 
data as shown in Figs 48 (Blattnig et al. 2005). Although 
mesons and their decay products play an important role in space 
radiation interaction with atmospheres where long interaction 
lengths allow decay into products which have small nuclear 
cross sections (Abe et al. 2003) and therefore buildup 
unhindered, the short interaction lengths in condensed materials 
as a spacecraft have only minor contributions to exposures from 
mesons (de Angelis et al. 2002). Most important to human 
space exploration will be their contribution on the mars surface 
which is now being addressed (Blattnig et al. 2005) for which 
the validation in Fig. 48 lends confidence in future HZETRN 
estimates of exposures on the Mars surface. 
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FUTURE ACTIVITY 

The early portion of this activity focused on the multiple 
charged ion shield code development including the 
atomic/nuclear database and transport procedures. As the early 
versions of these HZE transport codes matured (Wilson et al. 
1987a, 1987b) a separate code for proton/neutron transport was 
developed by a grant under the Historically Black College and 
University grant program at Hampton University. This code 
first needed a realistic nuclear database (Wilson et al. 1988a, 
1988c) for the development of computational procedures that 
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Fig. 48. Calculated muon spectra (top panel p + , 
bottom panel p") for three heights in Earth’s 
atmosphere compared with experimental results 
from BESS2001 (Abe et al. 2003). 
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could be seamlessly integrated with the HZE transport code using marching procedures (Wilson et al. 1991). With 
the improvement of the light ion transport procedures including the cluster knockout processes (Cucinotta 1993) the 
first generation of HZETRN was finalized (Wilson et al. 1995b) with the preliminary proton/neutron and light ion 
database. Among various improvements for the HZETRN, we unsuccessfully proposed improving the 
proton/neutron database as a critical issue in space shielding estimates but the proposal was rejected. Even in this 
new round of research and development, only transport code method development was funded with no support for 
nuclear database development. Since transport code development requires an adequate nuclear database some 


nuclear model investiments are made on a time available 
basis. Such an example is the pion production database 
developed by Steve Blattnig and John Norbury at the 
University of Wisconsin and the associated modified 
version of HZETRN (Blattnig et al. 2004) giving the 
results shown in Fig. 48. Under this new transport code 
development activity, the development of a fully 3-D 
computational procedure is being studied but with the 
preliminary nuclear models with minor modifications 
derived from “as-time-available” activity. A recent 
comparison of methods and databases is shown in Fig. 
49 including a forward/backward multigroup method 
(Clowdsley 2000), a Pn multigroup implementation 
(Heinbockel et al. 2003), an Sn coupled HZETRN 
approximation (Singleterry et al. 1998), and MCNPX 
solution for the neutron fluence generated by the 
February 23, 1956 solar particle event in aluminum. 
The forward/backward multigroup method (Fb) and 
polynomial expansion method (Pn) utilize the HZETRN 
database while the Sn and MCNPX utilize their own 
database. The fact that the similar results are found for 
the Fb and Pn methods indicates the most pressing 
problem is the nuclear database utilized in the HZETRN 
affecting the neutron transport. 


Neutron Fluence at depth of 50gm/cm 2 in 1 00 gm/cm 2 of Aluminum 
Shield exposed to the February 23, 1956 energic particle event. 



Fig. 49. Total neutron fluence at 50 g/cm2 depth of in a 
100 g/cm2 aluminum shield exposed to the fluence of the 
February 23, 1956 Solar particle event. 


SUMMARY 

In the present report, we have reviewed the status of the basic physics models and solution techniques for space 
radiation transport through materials and spacecraft of interest. This is followed by a review of the verification and 
validation (V&V) procedures used in code development. Since the V&V activity was applied during code as test to 
new sections of code, it is not a uniform work since differing stages of development were used. The final 
assessment of the marching version of the codes circa 1995 and later are mainly for the code as it stood after the 
end of the funding period for code development and so the code at that time has only seen changes in the neutron 
propagator. A satisfactory neutron propagator must await an improved nuclear data base development to be 
finalized. Even so the estimate of dosimetric quantities are judged to be uncertain by only 15 percent by workers at 
the Johnson Space Center and the basis of this assessment is partly presented in the validation work within this 
paper. Future work will emphasize both laboratory and space validation resulting in a flexible code useful not only 
to the design engineer but to the radiobiologist in their studies as well. 
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